function [U,L] = crout_decomposition(X)
U = zeros(size(X));
L = zeros(size(X));
for i = 1:length(X)
    L(1,i) = X(1,i);
end
for i = 1:length(X)
    U(i,1) = X(i,1)/X(1,1);
    U(i,i) = 1;
end
for j = 2:length(X)
    for i = 2:length(X)
        if j == i
            s = 0;
            for k = 1:j-1
                s = s + U(j,k)*L(k,j);
            end
            L(j,j) = X(j,j) - s;
        elseif j > i
            s = 0;
            for k = 1:i-1
                s = s + U(j,k)*L(k,i);
            end
            U(j,i) = (X(j,i) - s)/L(i,i);
        elseif j < i
            s = 0;
            for k = 1:j-1
                s = s + U(j,k)*L(k,i);
            end
            L(j,i) = X(j,i) - s;
        end
    end
end
disp(U)
disp(L)
end